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We use molecular dynamics simulations in 2d to study multi-component fluid in the limiting case 
where all the particles are different (APD). The particles are assumed to interact via Lennard-Jones 
(LJ) potentials, with identical size parameters but their pair interaction parameters are generated 
at random from a uniform or from a peaked distribution. We analyze both the global and the local 
properties of these systems at temperatures above the freezing transition and find that APD fluids 
relax into a non-random state characterized by clustering of particles according to the values of their 
pair interaction parameters (particle-identity ordering). 


Multi-component systems are often encountered in 
studies of colloidal fluids and dispersions, which are poly- 
disperse in particle size. They are also quite common in 
biology: a cell contains many thousands of different types 
of proteins and other macromolecules that differ in their 
size, shape and in their interactions. In this work we 
focus on the multiplicity of interactions between the con¬ 
stituents characteristic of biological systems (we make no 
attempt to model real proteins) and introduce a minimal 
physical model in which all particles are different (APD) 
in the sense that their interaction parameters are chosen 
at random from a given distribution. This model differs 
in several important ways from past theoretical studies 
of multi-component systems that addressed size polydis- 
persity of colloids and used a coarse-grained description 
(in terms of moments of the distribution of particle sizes) 
of the thermodynamic properties such as free energy and 
phase diagrams of these systems (see e.g., [IHS])- Con¬ 
versely, in our APD model all particles have identical 
size parameters and differ only in their attractive inter¬ 
actions with other particles. We use computer simula¬ 
tions to obtain insights about the temperature-dependent 
organization of APD fluids into non-random states in 
which the pair interaction parameters (PIP) of neigh¬ 
boring particles become strongly correlated. Because of 
the huge configurational space of APD systems (there are 
N\ non-equivalent particle permutations for each spatial 
configuration of N particles), we prefered to avoid apriori 
assumptions about equilibration and used molecular dy¬ 
namics (MD) rather than equilibrium Monte-Carlo meth¬ 
ods. 

The MD simulations were carried out in 2D under NVT 
ensemble {N = 2500 particles, density p* = iV/Area = 
0.6944) using the open source package LAMMPS [4]. All 
the particles have the same mass m and size cr which are 
set to unity, and the simulation box is periodic in both 
X- and Y- axes. The particles interact via Lennard-Jones 
(LJ) potential 

^L.iW = 4ezi [(cr/r)i2 _ : ( 0 - 1 ) 

which is cut and shifted to zero at a distance Vc = 2.5a 
(Clj(c > Tc) = 0). The pair interaction parameter 


gives the depth of the potential well and r is the sep¬ 
aration between a pair of particles i and j. All the 
physical quantities are expressed in LJ reduced units 
[^. The Langevin equations of motion of the particles 
are integrated with a timestep of St = 0.005Tj^j, where 
Tlj = a{m/eY/'^ = 1 is the LJ time (we take cr = to = 1 
and define e = 1 as the lower bound on e^). The friction 
coefficient is fixed at Q = l/rd = l/SOr^j (r^ is the vis¬ 
cous damping time). We equilibrate the system at high 
temperature and then bring the system down to the re¬ 
quired temperature and let it relax for about 5 x lO^r^j to 
steady-state in which the ensemble-averaged (potential) 
energy remains constant in time. All the data reported 
in this work were obtained in the stationary regime. 



FIG. 1: Heating curves of the potential energy per particle 
for the GM, U and IG (one-component) systems. The distri¬ 
bution of PIPs in the range 1 — 4 for GM and U systems is 
shown in inset (a), and the energy difference Aw between the 
IG-system and the APD-systems are shown in inset (b). 


Two different distributions P{eij) are considered: (1) 
Geometric mean (GM) distribution: a particle i is as¬ 
signed an interaction strength in the range 1 — 4 taken 
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randomly from a uniform distribution. We then adopt 
the Berthelot mixing ansatz and define Cij = 

The distribution P{cij) can be computed analytically us¬ 
ing the methods of ref. [7] and is shown in inset (a) in 
fig. flj It is peaked at e™'’ = 2.0 and has a mean of 
(ey) = 2.42. 

(2) Uniform (U) distribution: Instead of assigning each 
particle a value ey random ey values are assigned to all 
the possible pairs of particles, i.e., N{N — I)/2 numbers 
are drawn randomly from a uniform distribution in the 
range 1 — 4 (horizontal line in inset (a) in fig. [^. Unlike 
GM where the identity of a particle is determined by its 
value of ey in the U system a particle i is characterized 
by the set of — 1 PIPs with all other particles (j ^ i) 
and since the probability of generating 2 identical such 
sets is vanishingly small, we conclude that all particles in 
the system are different from each other (note that since 
their interaction parameters are taken from a uniform 
distribution, they are all identical in a statistical sense). 
As a reference we use a one-component (1C) system with 
Eij =2.5 (the average value in the interval 1 — 4) for all 
the pairs. 

In fig. we present the heating curves for the average 
potential energy per particle as a function of T for the two 
systems, GM and U (we used a heating-cooling rate of 
dT/dt « 10“®1/Tlj at which both systems exhibit weak 
hysteresis upon heating and cooling- not shown). The 
abrupt change of potential energy with temperature is a 
signature of a transition from a liquid to a solid phase, 
characterized by the appearance of positional and orien¬ 
tational (hexatic [8]) ordering. A detailed study of the 
low-temperature solid phase is beyond the scope of this 
work but preliminary results indicate that the particles 
form a hexagonal crystal (this crystal may be metastable 
with respect to PIP ordering). This concurs with the 
well-known observation that in order to prevent crystal¬ 
lization (which always occurs in simulations of 1C LJ 
systems), one has to resort to mixtures of particles of 
different sizes [5] (recall that our particles have identical 
size). The precise nature of the liquid-solid transition 
in 2D is still controversial even in 1C LJ systems and 
may be strongly affected by finite size effects [101. 
this work we concentrate on the statistical properties of 
APD fluids and, therefore, we focus on the temperature 
range T > T* where T* is defined as the temperature 
at which the radial distribution function decays over half 
the system size. This yields T* « 1.0 and 1.1, for the 
GM and the U systems respectively, the former being 
close to that of the 1C system. We verified that similar 
results are obtained if we identify the transition temper¬ 
ature with the peak of the specific heat, with the drop 
in mean inter-particle distance or with the increase in 
orientational ordering (not shown). Since the difference 
between the systems stems from different distributions of 
PIPs, we expect the energy vs. temperature curves to ap¬ 
proach each other in the large T limit (e^ <C ksT) where 


short-range (r < a) repulsion dominates, and to diverge 
from each other at lower temperatures {cij/ksT ~ 0(1)) 
where the non-universal attractive part of the potential 
energy plays an increasingly important role; both effects 
are observed in fig. Inspection of inset (b) in fig. 
shows that the GM system has a slightly higher energy 
than the 1C system, possibly because the mean of its 
P{€ij) distribution (2.42) is slightly lower than that of 
1C system (2.5). The larger slope of the potential en¬ 
ergy vs. temperature curve for the U system means that 
the specific heat Cv oc (f^)y and the entropy change 

AS = J^^{Cv/T)dT are both higher for the U system. 
If the energies of the U and the 1C systems were compa¬ 
rable, this would suggest that the liquid-solid transition 
should take place at lower temperature in the former sys¬ 
tem, in direct contradiction with our simulation results. 
This, however, is not the case and the upward shift of the 
transition temperature (1.1 for U vs. 0.97 for 1C) is the 
consequence of the large downward shift of the energy of 
the U system compared to the 1C observed in fig. As 
will be shown in the following, this shift reflects the PIP 
ordering of the U fluid due to the existence of a large 
number of states in which pairs of particles with high 
values of cluster together, thus reducing the mean en¬ 
ergy of the system. Although GM fluids undergo PIP 
ordering as well, its extent is more limited because of the 
constraints imposed by the peaked form of the distribu¬ 
tion function of the pair interaction parameters (see inset 
(a) in fig.0. 

We computed the probability distribution of the total 
energy and found that at high temperatures the distri¬ 
butions for all three systems are strictly Gaussian with 
width obeying the thermodynamic relation between the 
variance of energy fluctuations and the specific heat m- 
As the transition temperature is approached, deviations 
from Gaussian behavior are observed as the correlation 
length approaches the finite system size [T^, but while 
low energy fluctuations are enhanced and high energy 
ones are suppressed for the U fluids, the reverse takes 
place for GM and 1C fluids (see SI for plots of global 
energy distributions). This is another manifestation of 
higher PIP ordering in the U system. 

We proceed to take a closer look into the PIP ordering 
of APD fluids. In the GM fluid there is hierarchy in the 
spatial organization of particles according to Ci values, 
i.e., particles with larger values of tend to form clus¬ 
ters at intermediate T (around T*) while particles with 
small Ci are preferentially localized around vacancies in 
the dense phase (see fig. |^a)). This suggests that GM 
fluids relax to a locally-ordered state in which there are 
strong correlations between the interaction parameters of 
neighboring particles. In order to characterize the local 
statistical properties of both GM and U fluids (in the U 
system there is no single parameter that characterizes 
the i**' particle) we introduce the effective interaction pa- 
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FIG. 2: (a) Typical configuration of U and GM systems at 
T = 1.2 and T = 1.1 respectively, that correspond to 5 = 
{T — T*)/T* « 0.1. The particles are colored according to ef® 
values, (b) Probability density functions (PDF) of ef® for the 
APD systems. U system: Symbols are for the data obtained 
from tracking of three randomly chosen particles. Open and 
filled symbols represent two independent realizations of tij 
values. The red curve represents the ensemble average. GM 
system: Ensemble average data is shown. The inset shows 
the PDFs for particles with values in the range (1) 1-1.2, 
(2) 2.9-3.1 and (3) 3.8-4. 


shape that results from the constraints imposed by the 
non-uniform distribution P{eij). Further insight into the 
PIP ordering of the GM fluid can be gained by consid¬ 
ering particles with in a particular narrow interval of 
values and computing the distributions of e®® for three 
different such intervals. We find (see inset figj^b)) that 
the corresponding distributions are nearly Gaussian, with 
width increasing with . Since each of the peaks is much 
narrower than the cumulative distribution, the GM fluid 
behaves as a mixture of finite number (<C N) of com¬ 
ponents that undergoes microphase separation as T* is 
approached (see hgj^a)). 



FIG. 3: (a) T-dependence of the mean and variance of the 
distributions of ef® for the two APD systems, (b) The po¬ 
tential energy per particle as a function of the temperature, 
both scaled by (ei®)(r). 


rameter ef® of particle i, that depends not only on the 
intrinsic interaction parameters of the particle, but also 
on the identity of its nearest neighbors in a particular 
configuration of the system, 

rib 

ef = ^e„/n6, (0.2) 

t=i 

where the sum over j goes over all the nt surrounding par¬ 
ticles within a cut-off radius Cc = 1.7 that corresponds to 
the minimum between the first and the second peak of 
the radial distribution function. In figj^a) we show typi¬ 
cal snapshots of the two APD systems taken at the same 
normalized distance S = (T — t*)/T* k. 0.1 from their 
respective transition temperatures (the particles are col¬ 
ored according to the values of e®®) (see SI for movies). 
Inspection of this figure shows that the U system is more 
homogeneous and has a higher mean value of e®® than the 
GM system (for the 1C system the distribution is a delta 
function at e®® = 2.5). This observation is confirmed 
quantitatively in figsj^b) where we plotted the distri¬ 
butions of e®® for both systems. While the U distribu¬ 
tion is nearly Gaussian and is centered about e®® ~ 2.9, 
the GM distribution is much broader, with a complex 


We now turn to study the temperature dependence 
of the distributions of ef®. Inspection of figj^a) shows 
that at high temperatures, the means of the U and the 
GM distributions of e®® approach the means of the corre¬ 
sponding eij distributions, 2.5 and 2.42, respectively. As 
temperature is decreased the two distributions shift to 
higher values of e®® and, while the width of the U distri¬ 
bution decreases, that of the GM distribution increases. 
The upward shift of (e®®) is much more pronounced for 
the U than for the GM system, a signature of more ef¬ 
ficient PIP ordering in the former system. Preliminary 
results of our Monte-Carlo simulations of the crystaline 
phase of the U system indicate that as temperature is 
lowered below T*, the peak of the distribution contin¬ 
ues to shift to higher values of e®® and the distribution 
becomes asymmetric with a tail towards lower values of 
e®® (not shown). We believe that at yet lower temper¬ 
atures the distribution narrows down and approaches a 
(5-function as T —>■ 0 (in the limit N —?> cc). Since (e®®) 
can be thought of as new temperature-dependent energy 
scale, it is interesting to check whether the thermody¬ 
namics of the two multi-component systems can be re¬ 
duced to that of an effective one-component system [13] 
by rescaling all energies and temperatures by (e®®)(T). 
As shown in fig. [^b) this anzatz appears to work for 
both the U and the GM systems (in the latter case, some 
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qualitative differences with the 1C system remain upon 
rescaling). 

In order to gain additional insight into the nature of 
self-organization of APD fluids, we calculated the partial 
radial distribution function (PRDF) ^^(r) that gives the 
probability to find a particle of type j at a distance r 
from particle i, for GM and U fluids at their respective 
transition temperatures (see fig. |^. Here i stands for 
particles with ef® > (e|®) and j denotes either particles 
in the same (diagonal PRDF) or in the complementary 
range e®® < (e®®) (off-diagonal PRDF). Here (e®®) « 2.47 
(GM) and 2.97 (U). In the U system the only difference 
between the diagonal and the off-diagonal PRDFs is that 
the first peak is slightly higher for the former, in accord 
with our assertion that the PIP distribution in U fluids 
is spatially homogeneous. In the GM system, all the di¬ 
agonal peaks are higher than the off-diagonal ones, indi¬ 
cating that GM fluids relax into an inhomogeneous state 
in which regions enriched in particles with larger values 
of €i possess higher long-range (hexagonal) order. 



FIG. 4: Partial pair correlation function g.. (r) of GM and U 
systems (at their respective T*) obtained by gronping par¬ 
ticles according to their ef® values (see text). Diagonal and 
off-diagonal PRDFs are shown in red and green colors, re¬ 
spectively. 


We would like to stress that the distributions in 
figll^b), as well as all other statistical properties we con¬ 
sidered so far, do not depend on the particular realization 
of the chosen distribution, i.e., they are translation- 
ally ergodic m- In order to check whether APD fluids 
are ergodic in the usual sense of the word, we calculated 
the distribution of e®® by (a) sampling the instantaneous 
configuration of the system (ensemble average) and av¬ 
eraging over many instantaneous configurations in order 
to reduce the noise, and by (b) following the trajectories 
of several labeled particles over large time intervals and 


collecting the values of e®® along these trajectories. As 
shown in figj^b) both methods yield identical results for 
the U system (in the GM system trajectory sampling will 
yield different results for particles with different values of 
ei). 

To summarize, we have studied the steady-state prop¬ 
erties of systems of particles with random PIPs. We in¬ 
troduced an effective interaction parameter ef® that char¬ 
acterizes the degree of local PIP ordering and found that 
this parameter increases monotonically upon cooling as 
the freezing transition is approached (preliminary results 
of Monte-Carlo simulations of 2D APD crystals suggest 
that the increase of e®® saturates upon further lower¬ 
ing of temperature). The constraints that define this 
organized state confine the system to a subspace of the 
available configurational space but the dimensionality of 
this subspace appears to be sufficiently large to allow 
for its efficient exploration (for discussion of diffusion in 
rugged high-dimensional energy and htness landscapes 
see ref. m) and therefore, the macroscopic behavior of 
APD fluids is similar to that of normal fluids . The de¬ 
gree of PIP ordering depends on the distribution of PIPs 
and is much more pronounced in U than in GM systems. 
This is the consequence of the fact that the choice of 
PIPs is much more constrained in the GM than in the U 
system (we choose N random numbers for GM vs. 
random numbers for U), and due to the unbiased (flat) 
vs. peaked forms of the PIP distributions of U and GM 
systems, respectively. 

One may wonder whether APD systems obey standard 
thermodynamics since the fact that all particles are dis¬ 
tinguishable suggests that there is a non-extensive con¬ 
tribution to the entropy (InA^!) associated with all pos¬ 
sible permutations of the particles. However, in accord 
with general arguments about the resolution of the Gibbs 
paradox in classical systems of distinguishable colloidal 
particles misi, we did not observe any violations of stan¬ 
dard thermodynamic relations. Results on the dynamics 
of the systems, e.g., singe particle tracking (SPT) mea¬ 
surements which is more relevant to make connection 
with SPT experiments on complex (biological) systems 
(see, e.g., [nitts]) will be discussed in an upcoming pub¬ 
lication. Other distributions of pair interaction parame¬ 
ters and 3D APD systems will be studied in future work. 
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SUPPLEMENTARY MATERIAL 


Global energy fluctuations 


The distributions of the potential energy per particle u = U/N (with U the total potential energy, and N total number 
of particles) are obtained for temperatures above and close to T*. We rescale u to to get a distribution with zero 
mean and unit variance, i.e., we define a variable x = {u— {u))/au, with (u) the ensemble average and cr^ the variance 


and obtain the probability density function P{x) shown in fig SI 
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FIG. SI: PDF of X for the 1C and APD systems at three different temperatures indicated in the figure. The solid lines are the 
Gaussian PDF with zero-mean and unit variance. For temperatures approaching T* the PDF deviates from a Gaussian PDF. 



























